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Within Kirkwood theory, we study the translational diffusion coefficient of a single polymer chain 
in dilute solution, and focus on the small difference between the short-time Kirkwood value 
and the asymptotic long-time value D. We calculate this correction term by highly accurate large- 
scale Brownian Dynamics simulations, and show that it is in perfect agreement with the rigorous 
variational result D < D^ K ', and with Fixman's Green-Kubo formula, which is re-derived. This 
resolves the puzzle posed by earlier numerical results (Rey et al, Macromolecules 24, 4666 (1991)), 
which rather seemed to indicate D > D^ K ^ ; the older data are shown to have insufficient statistical 
accuracy to resolve this question. We then discuss the Green-Kubo integrand in some detail. This 
function behaves very differently for pre-averaged vs. fluctuating hydrodynamics, as shown for the 
initial value by analytical considerations corroborated by numerical results. We also present further 
numerical data on the chain's statics and dynamics. 



I. INTRODUCTION 

Diffusion of polymer chains in dilute solution has been 
a subject of substantial theoretical interest 1-11 . A partic- 
ularly important result is the Kirkwood formula for the 
translational diffusion coefficient of a molecule immersed 
in a solvent of temperature T and viscosity rj, 



D(K) = £o ksT /J_\ 
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Here fee denotes Boltzmann's constant, while N is the 
number of chain segments, each of which has a segmen- 
tal diffusion coefficient Dq = ksT/C,, where £ is the seg- 
mental friction coefficient, ( = 67rna, where a denotes 
the segment's Stokes radius. Rr is the hydrodynamic 
radius of the molecule, which is defined as a thermal av- 
erage over the inverse intramolecular distances between 
segments: 
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Equation 1 is not exact, but rather the result of a num- 
ber of approximations 11 . Firstly, the Brownian motion 
of the molecule is described by a Smoluchowski equation 
(Kirkwood's diffusion equation) in the space of the seg- 
ment coordinates f%. Secondly, the diffusion tensor Dij, 
whose off-diagonal elements describe the hydrodynamic 
interaction between segments i and j, is approximated 
by the Oseen formula, 
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where fij is the unit vector in direction of = r j — fj . 
Thirdly, the diffusion coefficient is evaluated in the short- 
time limit of the Smoluchowski equation. The diffusion 
coefficient in the asymptotic long-time limit, which we 
denote as D, is different from the short-time value D^ K \ 
due to the presence of intramolecular dynamic correla- 
tions, which are expected to decay on the time scale of 
the overall relaxation time of the chain's conformational 
degrees of freedom (the Zimm time Tz)- Thus, the be- 
havior on short time scales t <C Tz is governed by D^ K \ 
while on long time scales t 3> tz the diffusive motion 
takes place with the coefficient D. For a general diffu- 
sion tensor (not necessarily of the Oseen type) the same 
considerations apply. In this case the short-time diffu- 
sion coefficient is obtained by averaging over the trace: 
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in what follows, the term "Kirkwood formula" will always 
refer to Eq. 4. 

The difference between and D has caused con- 

siderable interest and some confusion (see below); it is 
the purpose of the present paper to contribute to the 
clarification of these questions. Both light scattering 
experiments 12-15 and first-principles Molecular Dynam- 
ics (MD) simulations 16,17 show that the difference must 
be quite small such that it is below resolution. In the case 
of MD, the main difficulty is the fact that the "short- 
time" limit of Kirkwood theory is of course not the ex- 
act t — > limit, but rather a time regime where the 
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local motions of the particles in their "cages" (or the mo- 
menta) have already relaxed completely, while the con- 
formational degrees of freedom of the chain must at the 
same time be considered as completely unrelaxcd. To 
reach this separation of time scales (i. e. ballistic / short- 
time diffusive / long-time diffusive) in an unambiguous 
fashion would require extremely long chains which arc 
not accessible to today's computer power. A related 
problem is the fact that the monomer diffusion coeffi- 
cient Do, which must be taken into account in order 
to evaluate correctly 17 , is a somewhat fuzzy and 

ill-defined object in MD. This is seen in a quite ob- 
vious fashion when looking at the velocity autocorrela- 
tion function of a single monomer, whose time integral 
is the chain diffusion constant D 18 : After a rather quick 
short-time decay, one gets a negative tail which slowly 
decays to zero on the time scale of the overall relax- 
ation of the chain as a whole 19 . This latter decay de- 
scribes the Rouse- or Zimm-like slowing down by the 
coupling to the other monomers, while the quick short 
time regime defines, in principle, the monomer diffusion 
coefficient D . However, there is no obvious and well- 
defined time which discriminates between short- and 
long-time regime. Depending on where this "cut" is be- 
ing put, one obtains slightly different values for Dq. This 
latter problem does not exist in a recently introduced hy- 
brid approach 20 , where the chain monomers are coupled 
dissipatively to a discretized hydrodynamic background, 
such that the monomer diffusion coefficient is an input 
parameter. However, the former problem (lack of sepa- 
ration of time scales) is still present. 

In order to cope with these difficulties, one therefore 
has to work strictly within the framework of Kirkwood 
theory (i. e. the Smoluchowski equation), such that Dq 
is just an input parameter, the ballistic regime does not 
exist, and the theory has a well-defined diffusive short 
time limit. Nevertheless, calculating the long-time dif- 
fusion coefficient D, and the smooth crossover between 
the short-time behavior governed by D^ K \ and the long- 
time regime governed by D, requires to actually solve the 
Smoluchowski equation. However, closed analytical solu- 
tions do not exist as soon as fluctuating hydrodynamic 
interactions and / or excluded volume interactions are 
introduced. For this reason, analytical progress on the 
problem has been limited; nevertheless, a few important 
results have been obtained. Within the preaveraged ap- 
proximation, Ottinger 21 derived a closed expression for D 
which clearly differs from D^ K \ Furthermore, it is possi- 
ble to establish rigorous variational bounds on transport 
coefficients 11,22,23 , from which one can show that D must 
be smaller than D^ K \ Moreover, Fixman 24 derived the 
Green-Kubo type relation 

D = D W - £>i (5) 

with 

^ 1 = 3A^iT dt ( /(0) ' 1W ) >0 (6) 



and 

A = 4 (?) 

ij 

where Fj is the force acting on segment j, and )Sy is the 

mobility tensor, /x^- = D^ / '(fc^T). 

However, for a more quantitative understanding one 
has to rely on numerical solutions of the Smoluchowski 
equation. The most common approach is Brownian Dy- 
namics (BD) simulations, which, for the case of fluctuat- 
ing hydrodynamic interactions, was pioneered by Ermak 
and McCammon 25 . BD simulations for polymer chains 
have been done by numerous workers before 24,26-36 . One 
of these studies 31 produced quite accurate data, which 
however seemed to indicate D > D^ K \ No satisfactory 
explanation of this deviation has been given so far. For 
this reason, we have looked into the problem again, and 
performed careful BD simulations of the same model with 
identical parameters. In this model, excluded volume in- 
teractions are represented by a relatively soft potential, 
while the hydrodynamic interactions are modeled via the 
Rotne-Prager-Yamakawa (RPY) tensor 22 ' 37 . 

As noted before 31 , this is a numerically demanding 
task, as such a calculation is looking for a rather small 
effect. One has, on the one hand, to watch out for dis- 
cretization errors due to the finite time step 31 , and, on 
the other hand, to make sure that the statistical accu- 
racy of the data suffices. In order to test the relation 
D = D^ — D\, all three quantities must be deter- 
mined independently. As our analysis shows (see below), 
it turns out that both and D\ can be calculated 

with moderate statistical effort, and Rcf. 31 has indeed 
obtained values which coincide with ours. The diffusion 
coefficient D itself, however, is subject to much larger 
fluctuations, such that its accurate determination needs 
substantially more computer power than what was avail- 
able to the authors of Ref. 31. We believe that this is 
the most likely explanation of the deviation, which is not 
confirmed by our much more accurate data (these are 
based on an observation time which is orders of mag- 
nitude larger than that of Ref. 31). Rather, we find a 
very nice agreement with the relation D = D^ K ^ — D\ 
with D\ > 0. We also re-derive Fixman's 24 Green-Kubo 
formula for D\, Eq. 6, by means of an alternative ap- 
proach based on the Mori-Zwanzig projection operator 
technique, and try to elucidate its relation to the stochas- 
tic analog of the velocity autocorrelation function. Fur- 
thermore, we attempt to discuss the scaling behavior of 
D\ via both scaling arguments and numerical data; here 
we have obtained some results which we view as some- 
what unexpected and counter-intuitive. However, these 
should be considered as preliminary and inconclusive. 
The computational demands of the calculations are so 
large that it is impossible to accurately study the dy- 
namics for chain lengths much beyond N ~ 10 2 , which is 



2 



clearly not long enough to justify any claim of asymptotic 
behavior. 

The outline of this article is as follows. In Sec. II 
we define the simulated model and describe the simu- 
lation procedure we applied. In Sec. Ill we derive the 
short-time diffusion coefficient and long-time diffusion 
coefficient from the mean square displacement of the cen- 
ter of mass, and compare the analytical predictions with 
the numerical results. Section IV presents our consider- 
ations concerning the scaling of the correction term D\ 
together with numerical data. Section V then discusses 
further data on the static and dynamic scaling properties, 
mainly to demonstrate the consistency of the results. In 
Sec. VI we conclude with some final remarks. 



where h is the time step, Fi the force on monomer i, 
Jtij = Dij/ '(ksT), and pi a random displacement with 
zero mean and variance-covariance matrix given by 

(p i ®p-' j } = 2D ij h. (10) 

We generated the random terms pi from a uniform 
distribution 38 , and used the procedure given in Ref. 25 
to satisfy Eq. 10. The term h^2- (d/drj) ■ on the 
right hand side of Eq. (9) was omitted, because the 
RPY tensor, due to the incompressibility of the solvent, 

is divergence-free, (d/drj) ■ = 0. 



II. MODEL AND SIMULATION METHOD 



D. Simulation Details 



A. Bead— Spring Model 

Following Rey et al. 31 , we have modeled the flexible 
polymer as a linear chain constituted by N units rang- 
ing from N = 6 to 200, each one connected with its first 
neighbors by means of harmonic springs. Therefore the 
bond lengths follow a Gaussian distribution whose vari- 
ance we denote by b 2 , i. e. b is the statistical segment 
length. Excluded volume forces only act between non- 
neighboring units; they are introduced by means of a 
potential of the form Aexp(—(3nj), where A and (3 are 
constant parameters. This potential is cut off at a dis- 
tance r r . 



B. Hydrodynamic Interaction 

The hydrodynamic interaction is introduced through 
the diffusion tensor proposed by Rotne and Prager and 
Yamakawa (RPY) 22,37 , 
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where a is the Stokes radius of the beads. This tensor is 
positive-definite for all chain configurations. 



C. Algorithm 

The polymer motion is governed by the stochastic dif- 
ferential equation that we solve through the first-order 
Ermak and McCammon algorithm 25 , 



\(t + h) = n(t) + ^2 V-a ■ Fj h + Ph 



(9) 



As in Ref. 31, we use a unit system where the three 
quantitities b, ksT, and £ = Qurja are set to unity. Hence 
the time unit is given by (b 2 /ksT. In this unit system, 
the other parameters are A = 75, /3 = 4, r c = 0.512, and 
a = 0.256. We used a reduced time step of h = 0.005, 
which is twice as small as that of Ref. 31. This choice 
was motivated by a comparison of end-to-end distance 
and gyration radius data obtained by simulations with 
various h with the corresponding data from Monte Carlo 
runs. This comparison showed that highly accurate re- 
sults require a rather conservative h value. 

For N < 50, data with good statistics were generated 
by the following procedure: (i) Equilibration by BD with- 
out hydrodynamic interactions, (ii) short additional equi- 
libration with hydrodynamic interactions turned on, and 
(iii) a long BD production run. Furthermore, we aver- 
aged over five independent such runs. For these chains, 
the resulting total simulation time divided by the longest 
relaxation time tz (which can be obtained from the time 
autocorrelation function of the end-to-end distance) is 
0.5 x 10 6 for N < 35, 0.4 x 10 6 for N = 40, 0.3 x 10 6 
for N = 45, and 0.2 x 10 6 for N = 50. All the statis- 
tical error bars were estimated by the blocking method 
developed by Flyvbjerg and Petersen 39 . 

From the results for these chains, it turned out that 
the dynamics on short time scales significantly below tz 
is of particular interest (see below). For this reason, we 
studied two longer chains (N = 100 and 200) by a some- 
what different procedure: By an efficient Monte Carlo 
procedure (the pivot algorithm 40 ' 41 combined with local 
moves) we generated a large sample of statistically inde- 
pendent conformations (128 conformations for N = 100, 
3700 conformations for N = 200), from which we started 
BD runs with hydrodynamic interactions of short to mod- 
erate length (in our time units: r = 2.5 x 10 4 for TV = 100, 
t = 100 for N — 200). Although this implies a quite good 
statistical accuracy, it is nevertheless significantly worse 
than for the shorter chains. Therefore only the most 
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important and interesting properties were evaluated for 
N = 100 and 200. 

For chain length N — 6, our program needed 60 sec. 
to run 10 5 BD steps on a 667 MHz PC. This number of 
course strongly increases with N (~ 500 sec. for TV = 50, 
~ 10 5 sec. for N = 200); the asymptotic iV 3 scaling of 
the algorithm was observed for (roughly) N > 30. The 
calculations ran for several months on a 16-machine PC 
cluster; the overall effort of the project in terms of single- 
processor time is estimated as roughly 6.7 years. 

III. TRANSLATIONAL DIFFUSION 
COEFFICIENT 

A. Theory 

In order to calculate the translational diffusion coeffi- 
cient in the short- and long-time regimes, we study the 
mean square discplacement of the center of mass, 

«cM = ^ r ; (ii) 

i 

From the Ermak-McCammon algorithm, Eq. 9, we find 
the updating rule for Rcm in one time step: 

AR(t) = R CM (t + h) - R CM {t) 

= ±(hA(t) + h^B(t)), (12) 

where we have introduced the abbreviations 

A = J2Vij-Fj (13) 

ij 

(this definition is identical to the notation of Fixman 24 ) 
and 

£ = /*- 1/2 X> (14) 

i 

Note that A and B are defined in such a way that they 
are independent of the time step h. As A and B, evalu- 
ated at the same time, are uncorrelatcd, and since (^B^ 
vanishes, we find 

= f^I>( B y} + (' ,2 >- < 15 > 

ij 

Thus the short-time diffusion coefficient is just given by 
the Kirkwood formula: 

DShort =^(M 2 ) 

= 3^£^)^ W - (16) 



For longer times (n time steps) we evaluate the mean 
square displacement as 

(R C M{nh) - R C m(0)) 2 ^ (17) 

= / AR(kh)\ )=J2 ( A R( kh ) ■ &R{lh)) ■ 
\ \fe=o Jim 

The matrix with elements I^AR{kh) ■ AR(lh)^ is obvi- 
ously symmetric. Since all elements with constant k — I 
are identical, for reasons of translational symmetry in 
time, we can simplify the previous expression as 

(^RcM{nh) - i? CM (0))^ (18) 

n-l 

= n (AR(0) 2 ^ + 2 J2( n - k ) ( A -R(°) ' &R(kh)^ . 
k=i 

This is quite analogous to the standard relation between 
mean square displacement and velocity autocorrelation 
function 18 . In the long-time limit n — > oo we thus ob- 
tain a diffusion coefficient which still depends on the time 
step: 

D ion g(h) = JL (^R CM (nh) - Rcm(0)) 2 ^ 

1 OO 

+ ^E( Ai? (°)- Ai? ( fc/i )); 

fe=i 

here we have assumed, as usual, that the correlation func- 
tion decays quickly to zero. 

We now make use of the fact that the stochastic terms 
B are uncorrelated at different times, and that there is a 
correlation between A and B only if A is evaluated at a 
later time as B (of course, a stochastic "kick" at a certein 
time will influence how the system evolves dynamically 
in the future). We thus obtain for k > 1 

(AR(0) ■ AR(kh)) = ^ (1(0) • A(kh)) 

+ ^(B(0).A(kh)), (20) 

yielding the relation 

D lon9 {h) = D (K) +D 1 {h)+D 2 {h) (21) 

with 

fc=i 
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and 
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In the continuum limit h — > 0, wc obviously have 
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(24) 



where the previous formulae tell us how the integral 
should be consistently discretized. Concerning D 2 , one 
obtains 



Do = 
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dt(B(0)-A(t) 



(25) 



At first glance, this looks as if this contribution would di- 
verge for h — > 0, but this is not the case. Rather, the cor- 
relation function (^B(0) ■ A(t)^ depends on the time step, 

and is, to leading order, proportional to h 1 / 2 , such that 
D 2 converges to a well-defined non-trivial value. This 

is demonstrated in Fig. 1, where ^-B(O) • A(t)^> h^ 1 ! 2 is 

plotted for N — 6 and various time steps h. The h 1 ^ 2 
dependence may be explained by linear response theory: 

As (A) vanishes for symmetry reasons, only that part 



of A will contribute to ^-B(O) • A{t)j which is actually 

the response to the "kick" at time zero. This "kick", 
however, has an infinitcsimally small amplitude of order 
h 1 / 2 . Therefore, linear response theory should be appli- 
cable, and the response in A should be proportional to 
h 1 / 2 as well. 

On a more formal level, we can write the Ermak- 
McCammon updating rule (Eq. 9) in the continuum limit 
as a Langevin equation with Ito interpretation: 
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Fj + fi 



(26) 



only conclusion is that D\ and D 2 must satisfy the rela- 
tion D 2 = —2D \. Unfortunately, we have not been able 
to derive this result directly; however, it is very nicely 
borne out by our numerical data (Table I). Furthermore, 

if (^A(0) ■ A(t)^ is positive for all times (as it is the case 

for our simulation data), it is immediately obvious that 
D must indeed be smaller than D^ K \ 

For the Mori-Zwanzig analysis of D, wc first notice 
that the Langevin equation corresponds to the Fokker- 
Planck equation (Kirkwood diffusion equation) 

J^p (r, t\r°, o) = -iCP (r, t\r°, o) , (29) 

where T is a shorthand notation for the set of all monomer 
coordinates, P(r,t|r°,0) is the transition probability 
density for the system going from T at time to T at 
time t, and — iC is the Fokker-Planck operator 



iC 



(30) 



where (3 = 1/(/cbT). The formal solution is P = 
exp(—iCt)S(T — r°), while the equilibrium distribution 
(i. e. the t — > oo solution) is 



P(r) 



exp(-0U) 
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(31) 



where U is the potential energy, such that Fi = —dU/dfi, 
and 
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for an arbitrary function f(T). Except for — iC, we 
also need iC\ which is the adjoint operator of — iC 
with respect to the standard scalar product (f\g) — 
J dTf(T)*g(r) (/* denoting the complex conjugate), 



where the noise term fi has zero mean, and variance 

[fi{t)®f j {t'))=2D i:j 5{t-t'), (27) 



i. e. /, corresponds to Pi/h. This in turn implies 
B = h 1 / 2 J2i fi, such that D 2 can also be written as 



Do 
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dt(fi(0)-A(t) 



(28) 



In seeming contrast to this result (D = D^ +Di +D 2 ), 
Fixman 24 rather obtained from linear response theory 
D = D( K ~) — D\, where D\ is defined precisely as in Eq. 
24. His result is however as valid as ours, and in what 
follows we will give an alternative derivation based on 
the Mori-Zwanzig projection operator formalism 42 . The 
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as well as —it, which is the adjoint operator of i& with 
respect to the natural scalar product 



(f\g) 



J dTp(T)f(Ty 9 (T). (34) 



From Eq. 32, and partial integration, one finds that — it 
coincides with iC) . 

We now use the standard memory equation as derived 
in Ref. 42. In Ref. 43 it was shown how to generalize this 
to the case of non-Hamiltonian dynamics; specifically it 
was shown there that the integral over the time correla- 
tion function of a slow variable S satisfies the relation 
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dt(s(oys(t)) = - (s\s) 2 x 



(35) 



(S\itf\S)+ / dt(S\idQexp(i£h)Qitf\S) 
Jo 



Here Q is the operator which projects onto the orthog- 
onal space of the slow variable, i. e. onto the space of 
variables which are (statically) uncorrelated with S. For 
the analysis of diffusion we study the variable 



S = cxp(«<f • R C m) 



(36) 



in the limit q — > such that q 1 is much larger than the 
polymer gyration radius. Therefore, 



(s(ors( t )) 



cxp 



iq-(Rc M (t)-Rc M (0)) ). (37) 



For times of order of the Zimm relaxation time, or 
smaller, this correlation function is very close to unity, 
due to the smallness of q. For times much larger, the 
motion of Rcm is just a Gaussian random walk with dif- 
fusion constant D, and hence 



{S(0yS(t)) = cxp \-!L ((R CM (t) - R CM (0)) 2 
= cxp(-Dq 2 t). 
Thus the left hand side of Eq. 35 is just 

1 



dt(S(0)*S(t)) 



Dq 
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(38) 



(39) 



As (S\S) = 1, we have 
D = -q- 2 (S\itf\S) 



dt(S\itfQexp(itft)Qitf\S) (40) 



in the limit q 
yields 



N 



= -jpY,4- 3 «-0:S+%q.AS 



iq 
N 



iq 



N' 



and, for q — > 0, 
<S|i£t|S> = 



f D (K) 



(41) 



(42) 



(the A term vanishes upon averaging, for symmetry rea- 
sons). This also means that in the limit q — > the vari- 
able iC^S becomes orthogonal to S, implying that in this 
limit we can ignore the operator Q in Eq. 40. Further- 
more, in the memory integral it is sufficient to just take 
the term linear in q for itfS — any higher order would 



not contribute to D in the limit q — > 0. In this order we 
find iC^S = iq ■ A/N. As — it = iC\ the memory term 
becomes 
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dt(S\itfQcxp(i£h)Qitf\S) 



dtl 1(0) -A(t)) =q 2 D 1 . 



(43) 



Combining these results, we obtain Fixman's 24 formula 
D = D<*)-Di. 



B. Numerical Results 

As we have seen in the previous subsection, the rela- 
tion 2Di + D 2 = should hold. Our data (see Table 
I) indeed confirm this prediction. Interestingly enough, 
we have found numerically that even the Green-Kubo 
integrands satisfy the corresponding relation 



A(t) =2(A(0)-A(t) 



B(0)-A(t)) = 0. (44) 



-1/2 



More precisely, we observed that for finite time step h 
there is a slight systematic deviation (A(i) ^ 0), which 
however tends to zero for h — > 0. Furthermore, we found 
that A(t) quickly decays to zero with increasing time t, 
and has both a positive (small t) and a negative (larger 
t) contribution, such that J Q dtA(t) is very small. Such 
discretization effects are the reason for our finding that 
the two functions 



and 



D(t) =£>W +Dx(t) +D 2 (t) 



D(t) = DM - £>!(*), 



0. Now, straightforward evaluation with 



Di(t) 
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D 2 (t) 
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dr (B(0) ■ A(t) 



(45) 



(46) 



(47) 



(48) 



are slightly different, in particular for short times. This 
is seen in Fig. 2, which shows D(t) and D(t) for N = 6, 
25, and 50, as a function of 1/t, with logarithmic ab- 
scissa. This figure also demonstrates that our data are 
well-converged and accurate enough to clearly discrimi- 
nate between short- and long-time regimes. 

Table I summarizes our results for the diffusion coef- 
ficient, where we list D<- K \ £>i and D 2 . The data con- 
firm the relation 2Di + D 2 = within our error bars up 
to N — 50. As it turned out that D 2 is much harder 
to sample than D\, we did not test the relation for the 
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longer chains N — 100 and 200, where our statistics is 
not sufficient. The data also show that the relative con- 
tribution of the correction term systematically increases 
with chain length (roughly 1% for TV = 6, roughly 3.5% 
for N = 200). 

Comparing our values for and D\ with those of 

Ref. 31, we see that they are compatible within error 
bars. However, Rey et al. 31 have obtained values for D 
which are larger than D^ K \ In view of this puzzle, we 
have done a test run for N — 6, where we increased the 
time step to their value h = 0.01, and decreased the ob- 
servation time to theirs (0.2 x 10 6 time steps). Figure 3 
shows our results for D\(t) and D 2 (t). One sees that the 
statistical accuracy is sufficient to obtain an acceptable 
value for D\ , but that it is by far not enough to estimate 
D 2 - We believe that this is the most likely explanation 
for the deviations observed in Ref. 31. 



IV. SCALING OF THE (AA) CORRELATION 
FUNCTION 

The systematic increase of the ratio D\/D^ K ^ with 
chain length, as seen from the data in Table I, raises 
the question if that ratio will saturate at a finite value, 
or keep on increasing, or maybe even tend to zero for 
N — > oo, after going through a maximum. We cannot 
give a conclusive answer to this question; however, we 
have found some interesting results concerning the issue. 
Rewriting the Green-Kubo formula for D\ as 

Dl= 3N*J dt ° A{tl m 

where Ca is the normalized A- A autocorrelation func- 
tion, 

C A (t) = j^(l(0).A(t)), (50) 

one sees that the N dependence is clear if it is known 
for (A 2 ) and for ta = f °° dt Ca (*) ■ As (A 2 ) is a static 
average, let us discuss it first. 

For fluctuating hydrodynamic interactions, we notice 
that Y]j Hij ■ Fj is nothing but the velocity flow field 
generated at position of monomer i, due to all the forces 
acting on the other monomers j. However, the system 
is in thermal equilibrium. Therefore, one should expect 
that this velocity is of order of a typical thermal velocity. 
Furthermore, in equilibrium the flow velocities at differ- 
ent volume elements are statistically uncorrclatcd. This 
picture suggests that A is essentially the sum of N statis- 
tically independent random variables, each of which does 
not depend on N. Therefore, the scaling 

(A 2 )oc7V (51) 

is expected from standard statistics, and this argument 
should be true independently of the details of the chain 
conformations. 



For preaveraged hydrodynamic interactions, however, 
this argument does not hold (the preaveraging prevents 
from "thermalizing"). Here we rather write 

W = EE(^)^( F M- (52) 

ij kl 

where summation over repeated Cartesian indices, de- 
noted by the Greek letters, is implied. Exploiting the 

isotropy of the (^P eg) F^ tensor (i. e. its proportionality 

to the unit tensor) one finds 

W = lEE(^)W(^^)' (53) 

ij kl 

This is simple to study for the case of a Gaussian 
chain, since then the (^F ■ F^ correlation is strictly short 

ranged. Indeed, for a random- walk chain, ^Fj ■ Fj^j must 

be zero if i and j are sufficiently far away from each other, 
since in that case one can choose a "pivot" monomer be- 
tween i and j, and rotate the "right" part of the chain 
around that monomer by a random angle, without chang- 
ing the statistical weight of the conformation. If j is on 
the rotated part, Fj is changed, while Fj is unchanged. 

Thus one shows (Pi ■ Pj} = - (Pi ■ P^ = 0. This 

argument holds whenever it is possible to find a pivot 
monomer, i. e. for \i — j\ > 2. Thus the only remain- 
ing correlations are those for i = j and \i — j\ = 1, in 
which case the correlation is obvious, due to the spring 
interaction with the neighboring monomer. In case of 
an excluded-volume chain we rather expect a power- 
law decay 44 , related to the probability of loops of length 
\i — j\. In what follows, we will therefore, for simplicity, 
focus on the Gaussian case. Noticing (/%) cx \i — jp 1 ^ 2 , 
we thus find 

pN pN pN 

(A 2 } cx / dx dy dz \x — y[' 1 ^ 2 \z — yP 1 ^ 2 
Jo Jo Jo 

(54) 

(the short range of (^F ■ F^ reduces the number of inte- 
grations from four to three). A trivial transformation to 
reduced variables x/N etc. then shows 

(A 2 ) oc iV 2 (55) 

for preaveraged hydrodynamics in the Gaussian case. 

We have tested these predictions numerically, and ex- 
ploited the fact that (^4 2 ) is a static average, and, as such, 
amenable to efficient Monte Carlo procedures. This is 
particularly true for the Gaussian case, where one simply 
generates a sample of chains. We were therefore able to 
study this case up to chains of length N — 0.8xl0 5 . How- 
ever, we restricted ourselves, for simplicity, to Oseen-like 
hydrodynamic interactions, where we studied both the 
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fluctuating and the preaveraged case. Apart from this, 
we also studied the behavior for our model (fluctuating 
hydrodynamics, excluded-volume chains) up to chains of 
length N = 10 4 . In this case, we used the full RPY 
interaction, and generated the conformations by a com- 
bination of the pivot algorithm 40,41 with local moves. For 
every chain length, 0.2 x 10 5 pivot moves, and 100 times 
as many local moves, were used. The results are pre- 
sented in Fig. 4; indeed reasonable agreement with our 
predictions is found. 

The scaling laws for (A 2 ) have an interesting implica- 
tion for the dynamics of A. Writing (A 2 ) oc N x where 
x = 1,2 for the discussed cases, and ta oc N v , we find 
from Eq. 49 



£>i oc N x+y - 



(56) 



On the other hand, it is well-established that is 
proportional to N'" where v is 1/2 for Gaussian chains, 
and 0.59 for excluded volume chains. It also strongly be- 
lieved that this is the asymptotic scaling law for D. This, 
however, implies that D\ must decay sufficiently quickly 
as a function of N — otherwise D\ would ultimately 
dominate and spoil the scaling of D. More precisely, one 
expects Di cx with (f)>v. Combined with the pre- 
vious consideration, this yields <j) — 2 — x — y > v or 
y < 2 — x — v, i. e. y < 1 — v for fluctuating hydrody- 
namics, and y < — v for preaveraged hydrodynamics of a 
Gaussian chain. This is a quite counter-intuitive result, 
since it implies that ta would increase only very weakly 
with chain length for fluctuating hydrodynamics, while 
it would even decrease for preaveraged hydrodynamics! 
Naively, one would rather expect that A, as a collective 
quantity, decays on the same time scale as the overall 
polymer conformations, i. e. ta oc tz oc N 3u (this latter 
relation is the standard Zimm scaling law 11 , and implies 
a rather sharp increase with N). We have thus found a 
"collective" quantity which apparently decays much more 
rapidly than the chain as a whole. We believe this is- 
sue deserves further attention; in particular, we think it 
would be very desirable to try to understand the under- 
lying physical mechanisms governing the relaxation of A 
somewhat better. 

Our numerical data from the BD simulation (i. c. 
for fluctuating hydrodynamics, and excluded-volume 
chains) can only give us very limited hints on the be- 
havior of ta as a function of N, since, due to the overall 
computational demand, we were not able to simulate the 
dynamics with sufficient accuracy for chains longer than 
N = 200. Our data for CU(i) are presented in Fig. 5. 
Apparently the correlation function has two distinct time 
regimes. In the short-time regime (t < to), the curves 
are practically superimposable. This is similar to the ob- 
servations made by Fixman 26 . In the long-time regime 
(t > t ), the correlation functions decay exponentially 
with a correlation time Trj, C^(t) oc exp(— t/Tjj). Figure 
6 shows our data for td- Indeed td increases with chain 
length; however, the observed behavior in our limited N 



window is anything but a power law. According to our 
previous considerations, the increase of td should not be 
stronger than N 1 ^". Indeed this condition seems to be 
satisfied in the regime of longer chains. 



V. FURTHER RESULTS 

A. Static Scaling Properties 

The radius of gyration and end-to-end distance are 
given by 



(57) 



and 



(Rl) = ((r N -nf). 



(58) 



The theoretical scaling for these static properties is 
(R 2 g )<x(R 2 e )oc(N-l) 2 \ 



(59) 



In good solvent, the scaling exponent has the theoreti- 
cal value of v « 0.588 from renormalization group cal- 
culations and Monte Carlo simulations 45 . The log-log 
fits of (i? 2 ) and (R 2 ,) vs. N — 1 yield the exponents 
2u = 1.187 ± 0.003 and 2v = 1.133 ± 0.006, respectively 
(see Fig. 7), which is similar to the results by Rey et 

r,l 3i 



Similarly, the static structure factor 

s w = ixi( exp ^-rij) 



N 



sin(fcrjj) 
kvij 



(60) 



which is measured in scattering experiments, obeys the 
scaling relation 



S(k) oc k- 1 ^ 



(61) 



in the regime R~ 1 <C k <C b^ 1 . By fitting a power law to 
our data we get the value (see Fig. 8) v = 0.575 ± 0.004. 

We have also obtained the first cumulant (or initial de- 
cay rate), Cl(k), of the dynamic structure factor S(k,t), 
defined as 



lim 



d ( S(k,t) 



t^o dt \S(k,0) J ' 



(62) 



Akcasu et al. 12 ' 46 have shown that fi(fc) can be written 
as 



Eij (k-Dij ■ kexp(ik-fij) 
n{k) = r '-. (63) 



E« (exp(ik ■ fij] 
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The orientational averaging in Eq. (63) is easily done for 
the RPY tensor 31 ' 47 . For the denominator one obtains 



CXp(ifc • fij) 



sin z 



6ni]a 



with z — krij . In the numerator we find for i = j 

(k ■ Dij ■ k) exp(ik ■ fj ) ) = 

For i ^ j one obtains instead 

(k • ■ k) exp(ik ■ fij) ' 
ksT , 2 



(64) 



(65) 




(66) 



in the case of large distances > 2a, while 

'^{k-Dij ■ k)exp(ik ■ fj)] 
k B T , 



6iri]a 



_^ 3 tij \ sin z 
16 a I z 



3 Tij ( oi is : sin : 



(67) 



for Tij < 2a. We should mention that there are some 
typographical errors both in Ref. 31 and Ref. 47. The 
right-hand term of Eq. (63) is therefore directly calcu- 
lated from the trajectories. In the k — > limit, Q(k) 
reflects exclusively the translational motion contribution 
to the chain dynamics. Therefore, the Kirkwood formula 
can be recovered from the first cumulant as 



D (K ^ = limO(fc)/fc 2 . 



(68) 



We have obtained D^ K ^ through Eq. (68) from the inter- 
cept of a fitting of Vt{k) /k 2 vs. k in the k — > limit shown 
in Fig. 9. These values are exactly the same as those 
obtained from Eq. (4), which constitutes a further verifi- 
cation of the consistency of our numerical method. From 
Fig. 9, a universal dependence of the type Q(k)/k 2 oc k 
is also obtained in the scaling regime i? _1 <C k <C b^ 1 . 



B. Dynamic Scaling Properties 

An approximately exponential behavior of the time 

correlation function (^R e (t) ■ R e (0)^, where R e is the 

end-to-end vector, is observed (see Fig. 10). We have 
extracted the relaxation times (Zimm times) tz corre- 
sponding to this behavior, tz is related to the orienta- 
tional diffusion of the end-to-end vector. A log-log fit 



of tz vs. N yields a slope of 1.71 ± 0.01 (Fig. 11), which 
is close to the theoretial value 3u with hydrodynamic in- 
teractions. 

The dynamic structure factor 



S(k,t) 



N ^ 



cxp 



ik-(fi(t)-fj(0)) 



(69) 



is predicted to exhibit scaling behavior 11 if both wave 
number and time are in the scaling regime, i. e. R~ 1 <C 
k ^ b^ 1 and To <C t <C Tz, where To is the microscopic 
time and tz is the Zimm time, the longest relaxation 
time of the chain. Fig. 12 gives a nice data collapse for 
the expected form 



S{k,t) 
S(k,0) 



f (k 2 t 2/3 ) 



(70) 



in log-linear represention. 

These scaling results demonstrate the internal consis- 
tency of our simulation, and in all cases agreement with 
the pertinent theories and experimental results. 



VI. SUMMARY 



The present study has shown that Brownian Dynam- 
ics simulations are able to attack the problem of trans- 
lational diffusion of polymer chains with hydrodynamic 
interaction and excluded volume. It has also highlighted 
the necessity of substantial statistical effort in order to 
obtain reliable data. While the standard picture of static 
and dynamic scaling is reproduced, as in previous stud- 
ies, the novel aspect is the calculation of the diffusion 
coefficient to sufficiently high accuracy, such that the dif- 
ference between the short-time Kirkwood value and the 
asymptotic long-time value could be resolved unambigu- 
ously. The numerical data are in perfect agreement with 
the theoretical predictions, both concerning the short- 
time value, and the crossover to the long-time value de- 
scribed by Fixman's Green-Kubo formula. It turns out 
that the long-time value is a few percent less than the 
short-time value. For Fixman's Green-Kubo integrand 
we find two remarkable results, namely that its initial 
value behaves very differently for preaveraged vs. fluctu- 
ating hydrodynamics, and that the correlation function, 
though describing a global property of the chain, must de- 
cay substantially faster than the conformations, in order 
to avoid a violation of dynamic scaling. Our numerical 
data are in reasonable agreement with these considera- 
tions, but not fully conclusive since only short chains were 
accessible. More work on this issue, in particular aimed 
at a better physical understanding, is clearly desirable. 
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A T 

N 






7-1 




7-1 






2£>i + D 2 






6 


r\ or ,( ,\ 1 -i r\- 

0.3544 ± 10 


-4 


0.00408 ± 10 


-5 


—0.0081 ± 


10" 


-4 


0.0000 ± 10 


-4 


A A1 1 O I O 1 f \ — 4 

0.0113 ± 3 x 10 


8 


0.3011 ± 10 


-4 


0.00441 ± 10 


-5 


—0.0087 ± 


10" 


-4 


0.0001 ± 10 


-4 


0.0143 ± 3 x 10 


11 


0.2513 ± 10 


-4 


0.00450 ± 10 


-5 


—0.0089 ± 


10" 


-4 


0.0001 ± 10 


-4 


0.0175 ± 4 x 10 


15 


0.2106 ± 10 


-4 


r\ r\ / \ a Af\ i -i r\- 

0.00443 ± 10 




—0.0088 ± 


10" 


-4 


0.0000 ± 10 


-4 


0.0214 ± 5 x 10 


20 


0.1788 ± 10 


-4 


0.00443 ± 10 


-5 


—0.0084 ± 


10" 


-4 


0.0001 ± 10 


-4 


0.0234 ± 6 x 10 


25 


0.1573 ± 10 


-4 


r\ r\f\ A C\C\ 1 -1 r\- 

0.00422 ± 10 


-5 


A AA^7A 1 

—0.0079 ± 


10" 


-4 


0.0000 ± 10 


-4 


0.0249 ± 6 x 10 


30 


0.1416 ± 10 


-4 


0.00399 ± 10 


-5 


—0.0076 ± 


10" 


-4 


r\ aaai i i r\ — 

0.0001 ± 10 


-4 


0.0267 ± 7 x 10 


35 


0.1298 ± 10" 


-4 


0.00381 ± 10 


-5 


-0.0072 ± 


10" 


-4 


0.0000 ± 10" 


-4 


0.0275 ± 8 x 10~ 4 


40 


0.1201 ± 10" 


-4 


0.00363 ± 10 


-5 


-0.0069 ± 


10" 


-4 


0.0000 ± 10" 


-4 


0.0287 ±8 x 10~ 4 


45 


0.1123 ±10" 


-4 


0.00345 ± 10 


-5 


-0.0066 ± 


10" 


-4 


0.0000 ± 10" 


-4 


0.0292 ± 9 x 10~ 4 


50 


0.1055 ±10" 


-4 


0.00332 ± 10 


-5 


-0.0063 ± 


10" 


-4 


0.0000 ± 10" 


-4 


0.0298 ± 9 x 10" 4 


100 


0.0718 ± 10" 


-4 


0.0022 ± 10" 


4 












0.031 ± 1 x 10~ 3 


200 


0.0428 ± 10" 


-4 


0.0015 ± 10" 


4 












0.035 ± 1 x 10~ 3 



TABLE I. The diffusion coefficients D {K) , D lt D 2 , as well as 2D l + D 2 , and |£>i + D 2 \ /D {K) , as defined in the text, for 
different chain lengths N. Note that D 2 was not sampled for N = 100 and 200, for reasons of poor statistics, and that hence 
for these chains we have assumed \Di + D 2 \ = D\. 
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10" 1 10° 10 1 

t 

FIG. 1. The correlation function ^A(t) ■ B(0)}, scaled with h~ 1/2 , for N = 6 and various time steps h = 0.0005, 0.001, 0.002, 
0.005 and 0.01, respectively. 



12 
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0.155 ■ 
0.154 ■ 
0.153 




FIG. 2. Time-dependent diffusion coefficients D(t) = D (K) + Di(t) + D 2 (t) (dotted lines) and D(t) = D (K) - Di(t) (full 
lines; see also Eq. 45 etc.), for = 6 (upper graph), = 25 (middle graph) and N = 50 (lower graph). The dashed lines 
indicate the corresponding Kirkwood values D^ K >. 
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FIG. 3. Time-dependent diffusion coefficients -Di(i) and D2(t) for N = 6 for different overall observation times: (i) Solid 
line: D\(t) (high resolution), (ii) Solid line with triangles: £>2(£) (high resolution), (iii) Solid line with circles: Di(t) (low 
resolution), (iv) Dashed line: Z>2(i) (low resolution). The parameters of the low-resolution calculation are adapted to those of 
Rey et a/. 31 , as explained in the text. 
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FIG. 4. Static average (A 2 ) as a function of chain length N for different cases: (i) Gaussian chain with preaveraged 
hydrodynamics (circles); (ii) Gaussian chain with fluctuating hydrodynamics (triangles); (iii) excluded-volume chain with 
fluctuating hydrodynamics (squares). The slopes of the solid lines indicate the power laws N 2 and N 1 . 
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FIG. 5. Normalized autocorrelation function {A(t)A(0)J / \A 2 J, for N = 6, 8, 11, 15, 20, 25, 30, 40, 50, 100, 200. The chain 
length increases from left to right. 
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10 100 

N 

FIG. 6. Chain length scaling for td- The slope of the solid line indicates the power law iV 1 
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AM 

FIG. 7. Chain length scaling for (R 2 g ) and (R 2 e ). The error bars are smaller than the symbol size. 
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FIG. 12. Log-linear scaling plot of the decay of the dynamic structure factor for N = 15, 20, 25, 30, 35, 40, 45, 50 for Zimm 
scaling: S(k,t)/S(k,0) vs. k 2 t 2 ^ . We restricted the wave numbers to the values k — 1.0, 1.5, 2.0, and the time regime to 
5 < t < 20. 
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